home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / raytrace / rayshade / src / blob.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-28  |  18.5 KB  |  711 lines
  1. /*
  2.  * blob.c
  3.  *
  4.  * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: blob.c,v 4.0 91/07/17 14:36:07 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    blob.c,v $
  19.  * Revision 4.0  91/07/17  14:36:07  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "blob.h"
  25.  
  26. static Methods *iBlobMethods = NULL;
  27. static char blobName[] = "blob";
  28.  
  29. unsigned long BlobTests, BlobHits;
  30.  
  31. /*
  32.  * Blob/Metaball Description
  33.  *
  34.  * In this implementation a Blob is made up of a threshold T, and a 
  35.  * group of 1 or more Metaballs.  Each metaball 'i'  is defined by 
  36.  * its position (xi,yi,zi), its radius ri, and its strength si.
  37.  * Around each Metaball is a density function di(Ri) defined by:
  38.  *
  39.  *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
  40.  *         di(Ri) =  0                                    ri < Ri 
  41.  *
  42.  * where
  43.  *       c4i  = si / (ri^4)
  44.  *       c2i  = -(2 * si) / (ri^2)
  45.  *       c0i  = si
  46.  *       Ri^2 is the distance from a point (x,y,z) to the center of
  47.  *            the Metaball - (xi,yi,zi).
  48.  *
  49.  * The density function looks sort of like a W (Float-u) with the 
  50.  * center hump crossing the d-axis at  si  and the two humps on the
  51.  * side being tangent to the R-axis at  -ri  and  +ri. Only the 
  52.  * range  [0,ri] is being used. 
  53.  * I chose this function so that derivative = 0  at the points  0 and +ri. 
  54.  * This keeps the surface smooth when the densities are added. I also 
  55.  * wanted no  Ri^3  and  Ri  terms, since the equations are easier to 
  56.  * solve without them. This led to the symmetry about the d-axis and
  57.  * the derivative being equal to zero at -ri as well.
  58.  *
  59.  * The surface of the Blob is defined by:
  60.  *
  61.  *
  62.  *                  N
  63.  *                 ---  
  64.  *      F(x,y,z) = \    di(Ri)  = T
  65.  *                 /
  66.  *                 ---
  67.  *                 i=0
  68.  *
  69.  * where
  70.  *
  71.  *     di(Ri)    is   given above
  72.  *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
  73.  *      N        is   number of Metaballs in Blob
  74.  *      T        is   the threshold
  75.  *    (xi,yi,zi) is   the center of Metaball i
  76.  *
  77.  */
  78.  
  79. /*****************************************************************************
  80.  * Create & return reference to a metaball.
  81.  */
  82. Blob *
  83. BlobCreate(T, mlist, npoints)
  84. Float T;
  85. MetaList *mlist;
  86. int npoints;
  87. {
  88.     Blob *blob;
  89.     int i;
  90.     MetaList *cur;
  91.  
  92. /* 
  93.  * There has to be at least one metaball in the blob.
  94.  * Note: if there's only one metaball, the blob 
  95.  * will be a sphere.
  96.  */
  97.     if (npoints < 1)
  98.     {
  99.         RLerror(RL_WARN, "blob field not correct.\n");
  100.         return (Blob *)NULL;
  101.     }
  102.  
  103. /*  
  104.  * Initialize primitive and Geom structures
  105.  */
  106.     blob = (Blob *)RayMalloc(sizeof(Blob));
  107.     blob->T = T;
  108.     blob->list=(MetaVector *)
  109.         RayMalloc( (unsigned)(npoints*sizeof(MetaVector)) );
  110.     blob->num = npoints;
  111.  
  112. /*
  113.  * Initialize Metaball list
  114.  */
  115.     for(i=0;i<npoints;i++)
  116.     {
  117.         cur = mlist;
  118.         if ( (cur->mvec.c0 < EPSILON) || (cur->mvec.rs < EPSILON) )
  119.         {
  120.             RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n");
  121.             return (Blob *)NULL;
  122.         }
  123.         /* store radius squared */
  124.         blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
  125.         /* Calculate and store coefficients for each metaball */
  126.         blob->list[i].c0 = cur->mvec.c0;
  127.         blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
  128.         blob->list[i].c4 = cur->mvec.c0 /
  129.                 (blob->list[i].rs * blob->list[i].rs);
  130.         blob->list[i].x = cur->mvec.x;
  131.         blob->list[i].y = cur->mvec.y;
  132.         blob->list[i].z = cur->mvec.z;
  133.         mlist = mlist->next;
  134.         free((voidstar)cur);
  135.     }
  136.     /* 
  137.      * Allocate room for Intersection Structures and
  138.      * Allocate room for an array of pointers to point to
  139.      * the Intersection Structures.
  140.      */
  141.     blob->ilist = (MetaInt *)RayMalloc( 2 * blob->num * sizeof(MetaInt));
  142.     blob->iarr = (MetaInt **)RayMalloc( 2 * blob->num * sizeof(MetaInt *));
  143.     return blob;
  144. }
  145.  
  146. Methods *
  147. BlobMethods()
  148. {
  149.     if (iBlobMethods == (Methods *)NULL) {
  150.         iBlobMethods = MethodsCreate();
  151.         iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
  152.         iBlobMethods->methods = BlobMethods;
  153.         iBlobMethods->name = BlobName;
  154.         iBlobMethods->intersect = BlobIntersect;
  155.         iBlobMethods->normal = BlobNormal;
  156.         iBlobMethods->bounds = BlobBounds;
  157.         iBlobMethods->stats = BlobStats;
  158.         iBlobMethods->checkbounds = TRUE;
  159.         iBlobMethods->closed = TRUE;
  160.     }
  161.     return iBlobMethods;
  162. }
  163.  
  164. /*****************************************************************************
  165.  * Function used by qsort() when sorting the Ray/Blob intersection list
  166.  */
  167. int
  168. MetaCompare(A,B)
  169. char *A,*B;
  170. {
  171.     MetaInt **AA,**BB;
  172.  
  173.     AA = (MetaInt **) A;
  174.     BB = (MetaInt **) B;
  175.     if (AA[0]->bound == BB[0]->bound) return(0);
  176.     if (AA[0]->bound <  BB[0]->bound) return(-1);
  177.     return(1);  /* AA[0]->bound is > BB[0]->bound */
  178. }
  179.  
  180. /*****************************************************************************
  181.  * Ray/metaball intersection test.
  182.  */
  183. int
  184. BlobIntersect(blob, ray, mindist, maxdist)
  185. Blob *blob;
  186. Ray *ray;
  187. Float mindist, *maxdist;
  188. {
  189.     double c[5], s[4];
  190.     Float dist;
  191.     MetaInt *ilist,**iarr;
  192.     register int i,j,inum;
  193.     extern void qsort();
  194.  
  195.     BlobTests++;
  196.  
  197.     ilist = blob->ilist;
  198.     iarr = blob->iarr;
  199.  
  200. /*
  201.  * The first step in calculating the Ray/Blob intersection is to 
  202.  * divide the Ray into intervals such that only a fixed set of 
  203.  * Metaballs contribute to the density function along that interval. 
  204.  * This is done by finding the set of intersections between the Ray 
  205.  * and each Metaball's Sphere/Region of influence, which has a 
  206.  * radius  ri  and is centered at (xi,yi,zi).
  207.  * Intersection information is kept track of in the MetaInt 
  208.  * structure and consists of:
  209.  *
  210.  *   type    indicates whether this intersection is the start(R_START)
  211.  *           of a Region or the end(R_END) of one. 
  212.  *   pnt     the Metaball of this intersection
  213.  *   bound   the distance from Ray origin to this intersection
  214.  *
  215.  * This list is then sorted by  bound  and used later to find the Ray/Blob 
  216.  * intersection.
  217.  */
  218.  
  219.     inum = 0;
  220.     for(i=0; i < blob->num; i++)
  221.     {
  222.         register Float xadj, yadj, zadj;
  223.         register Float b, t, rs;
  224.         register Float dmin,dmax;
  225.  
  226.         rs   = blob->list[i].rs;
  227.         xadj = blob->list[i].x - ray->pos.x;
  228.         yadj = blob->list[i].y - ray->pos.y;
  229.         zadj = blob->list[i].z - ray->pos.z;
  230.  
  231.     /*
  232.      * Ray/Sphere of Influence intersection
  233.      */
  234.         b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
  235.         t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
  236.  
  237.     /* 
  238.      * don't except imaginary or single roots. A single root is a ray
  239.      * tangent to the Metaball's Sphere/Region. The Metaball's
  240.      * contribution to the overall density function at this point is
  241.      * zero anyway.
  242.      */
  243.         if (t > 0.0) /* we have two roots */
  244.         {
  245.             t = sqrt(t);
  246.             dmin = b - t;
  247.     /* 
  248.      * only interested in stuff in front of ray origin 
  249.      */
  250.             if (dmin < mindist) dmin = mindist;
  251.             dmax = b + t;
  252.             if (dmax > dmin) /* we have a valid Region */
  253.             {
  254.         /*
  255.              * Initialize min/start and max/end Intersections Structures
  256.              * for this Metaball
  257.              */
  258.                 ilist[inum].type = R_START;
  259.                 ilist[inum].pnt = i;
  260.                 ilist[inum].bound = dmin;
  261.                 for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  262.                 iarr[inum] = &(ilist[inum]);
  263.                 inum++;
  264.  
  265.                 ilist[inum].type = R_END;
  266.                 ilist[inum].pnt = i;
  267.                 ilist[inum].bound = dmax;
  268.                 for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  269.                 iarr[inum] = &(ilist[inum]);
  270.                 inum++;
  271.             } /* End of valid Region */
  272.         } /* End of two roots */
  273.     } /* End of loop through metaballs */
  274.  
  275.     /*
  276.      * If there are no Ray/Metaball intersections there will 
  277.      * not be a Ray/Blob intersection. Exit now.
  278.      */
  279.     if (inum == 0)
  280.     {
  281.         return FALSE;
  282.     }
  283.  
  284.     /* 
  285.      * Sort Intersection list. No sense using qsort if there's only
  286.      * two intersections.
  287.      *
  288.      * Note: we actually aren't sorting the Intersection structures, but
  289.      * an array of pointers to the Intersection structures. 
  290.      * This is faster than sorting the Intersection structures themselves.
  291.      */
  292.     if (inum==2)
  293.     {
  294.         MetaInt *t;
  295.         if (ilist[0].bound > ilist[1].bound)
  296.         {
  297.             t = iarr[0];
  298.             iarr[0] = iarr[1];
  299.             iarr[1] = t;
  300.         }
  301.     }
  302.     else qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *),
  303.             MetaCompare);
  304.  
  305. /*
  306. * Finding the Ray/Blob Intersection
  307. * 
  308. * The non-zero part of the density function for each Metaball is
  309. * 
  310. *   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i 
  311. *
  312. * To find find the Ray/Blob intersection for one metaball
  313. * substitute for distance   
  314. *
  315. *     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
  316. * 
  317. * and then substitute the Ray equations:
  318. *   
  319. *     x  = x0 + x1 * t
  320. *     y  = y0 + y1 * t
  321. *     z  = z0 + z1 * t
  322. *
  323. * to get a big mess :^). Actually, it's a Quartic in t and it's fully 
  324. * listed farther down. Here's a short version:
  325. *
  326. *   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
  327. *
  328. * Don't forget that the Blob is defined by the density being equal to 
  329. * the threshold T.
  330. * To calculate the intersection of a Ray and two or more Metaballs, 
  331. * the coefficients are calculated for each Metaball and then added 
  332. * together. We can do this since we're working with polynomials.
  333. * The points of intersection are the roots of the resultant equation.
  334. *
  335. * The algorithm loops through the intersection list, calculating
  336. * the coefficients if an intersection is the start of a Region and 
  337. * adding them to all intersections in that region.
  338. * When it detects a valid interval, it calculates the 
  339. * roots from the starting intersection's coefficients and if any of 
  340. * the roots are in the interval, the smallest one is returned.
  341. *
  342. */
  343.  
  344.     {
  345.         register Float *tmpc;
  346.         MetaInt *strt,*tmp;
  347.         register int istrt,itmp;
  348.         register int num,exitflag,inside;
  349.  
  350.     /*
  351.      * Start the algorithm outside the first interval
  352.      */
  353.         inside = 0;
  354.         istrt = 0; 
  355.         strt = iarr[istrt];
  356.         if (strt->type != R_START)
  357.             RLerror(RL_WARN,"MetaInt sanity check FAILED!\n");
  358.  
  359.         /*
  360.          * Loop through intersection. If a root is found the code
  361.          * will return at that point.
  362.          */
  363.         while( istrt < inum )
  364.         {
  365.             /* 
  366.              * Check for multiple intersections  at the same point.
  367.              * This is also where the coefficients are calculated
  368.              * and spread throughout that Metaball's sphere of
  369.              * influence.
  370.              */
  371.             do
  372.             {
  373.                 /* find out which metaball */
  374.                 i = strt->pnt;
  375.                 /* only at starting boundaries do this */
  376.                 if (strt->type == R_START)
  377.                 {
  378.                     register MetaVector *ml;
  379.                     register Float a1,a0;
  380.                     register Float xd,yd,zd;
  381.                     register Float c4,c2,c0;
  382.  
  383.                     /* we're inside */
  384.                     inside++;
  385.  
  386.     /*  As promised, the full equations
  387.      *
  388.      *   c[4] = c4*a2*a2; 
  389.      *   c[3] = 4.0*c4*a1*a2;
  390.      *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
  391.      *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
  392.      *   c[0] = c4*a0*a0 + c2*a0 + c0;
  393.      *
  394.      * where
  395.      *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray 
  396.      *                                           is normalized
  397.      *        a1 = (xd*x1 + yd*y1 + zd*z1)
  398.      *        a0 = (xd*xd + yd*yd + zd*zd)
  399.      *        xd = (x0 - xi)
  400.      *        yd = (y0 - yi)
  401.      *        zd = (z0 - zi)
  402.      *        (xi,yi,zi) is center of Metaball
  403.      *        (x0,y0,z0) is Ray origin
  404.      *        (x1,y1,z1) is normalized Ray direction
  405.      *        c4,c2,c0   are the coefficients for the
  406.      *                       Metaball's density function
  407.      *
  408.      */
  409.  
  410.                     /*
  411.                      * Some compilers would recalculate
  412.                      * this each time its used.
  413.                      */
  414.                     ml = &(blob->list[i]);
  415.  
  416.                     xd = ray->pos.x - ml->x;
  417.                     yd = ray->pos.y - ml->y;
  418.                     zd = ray->pos.z - ml->z;
  419.                     a1 = (xd * ray->dir.x + yd * ray->dir.y
  420.                         + zd * ray->dir.z);
  421.                     a0 = (xd * xd + yd * yd + zd * zd);
  422.  
  423.                     c0 = ml->c0;
  424.                     c2 = ml->c2;
  425.                     c4 = ml->c4;
  426.  
  427.                     c[4] = c4;
  428.                     c[3] = 4.0*c4*a1;
  429.                     c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
  430.                     c[1] = 2.0*a1*(2.0*c4*a0 + c2);
  431.                     c[0] = c4*a0*a0 + c2*a0 + c0;
  432.  
  433.         /* 
  434.          * Since we've gone through the trouble of calculating the 
  435.          * coefficients, put them where they'll be used.
  436.          * Starting at the current intersection(It's also the start of
  437.          * a region), add the coefficients to each intersection until
  438.          * we reach the end of this region.
  439.          */
  440.                     itmp = istrt; 
  441.                     tmp = strt;
  442.                     while( (tmp->pnt != i) ||
  443.                            (tmp->type != R_END) )
  444.                     {
  445.                         tmpc = tmp->c;
  446.                         for(j=0;j<5;j++)
  447.                             *tmpc++ += c[j];
  448.                         itmp++; 
  449.                         tmp = iarr[itmp];
  450.                     }
  451.  
  452.                 } /* End of start of a Region */ 
  453.                 /* 
  454.                  * Since the intersection wasn't the start
  455.                  * of a region, it must the end of one.
  456.                  */
  457.                 else inside--;
  458.  
  459.         /*
  460.          * If we are inside a region(or regions) and the next
  461.          * intersection is not at the same place, then we have
  462.          * the start of an interval. Set the exitflag.
  463.          */
  464.                 if ((inside > 0)
  465.                     && (strt->bound != iarr[istrt+1]->bound) )
  466.                     exitflag = 1;
  467.                 else 
  468.                 /* move to next intersection */
  469.                 {
  470.                     istrt++; 
  471.                     strt = iarr[istrt];
  472.                     exitflag = 0;
  473.                 }
  474.         /* 
  475.          * Check to see if we're at the end. If so then exit with
  476.          * no intersection found.
  477.          */
  478.                 if (istrt >= inum)
  479.                 {
  480.                     return FALSE;
  481.                 }
  482.             } while(!exitflag);
  483.                 /* End of Search for valid interval */
  484.  
  485.         /*
  486.          * Find Roots along this interval
  487.          */
  488.  
  489.             /* get coefficients from Intersection structure */
  490.             tmpc = strt->c;
  491.             for(j=0;j<5;j++) c[j] = *tmpc++;
  492.  
  493.             /* Don't forget to put in threshold */
  494.             c[0] -= blob->T;
  495.  
  496.             /* Use Graphic Gem's Quartic Root Finder */
  497.             num = SolveQuartic(c,s);
  498.  
  499.         /*
  500.          * If there are roots, search for roots within the interval.
  501.          */
  502.             dist = 0.0;
  503.             if (num>0)
  504.             {
  505.                 for(j=0;j<num;j++)
  506.                 {
  507.         /* 
  508.          * Not sure if EPSILON is truly needed, but it might cause
  509.          * small cracks between intervals in some cases. In any case 
  510.          * I don't believe it hurts.
  511.          */
  512.                     if ((s[j] >= (strt->bound - EPSILON))
  513.                         && (s[j] <= (iarr[istrt+1]->bound +
  514.                             EPSILON) ) )
  515.                     {
  516.                         if (dist == 0.0)
  517.                         /* first valid root */
  518.                             dist = s[j];
  519.                         else if (s[j] < dist)
  520.                          /* else only if smaller */
  521.                             dist = s[j];
  522.                     }
  523.                 }
  524.                 /*
  525.                  * Found a valid root 
  526.                  */
  527.                 if (dist > mindist && dist < *maxdist)
  528.                 {
  529.                     *maxdist = dist;
  530.                     BlobHits++;
  531.                     return TRUE;
  532.                     /* Yeah! Return valid root */
  533.                 }
  534.             }
  535.  
  536.         /* 
  537.          * Move to next intersection
  538.          */
  539.             istrt++;
  540.             strt = iarr[istrt];
  541.  
  542.         } /* End of loop through Intersection List */
  543.     } /* End of Solving for Ray/Blob Intersection */
  544.  
  545.     /* 
  546.      * return negative
  547.      */
  548.     return FALSE;
  549. }
  550.  
  551.  
  552. /***********************************************
  553.  * Find the Normal of a Blob at a given point
  554.  *
  555.  * The normal of a surface at a point (x0,y0,z0) is the gradient of that
  556.  * surface at that point. The gradient is the vector 
  557.  *
  558.  *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
  559.  *
  560.  * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
  561.  * to x,y and z respectively. Since the surface of a Blob is made up
  562.  * of the sum of one or more polynomials, the normal of a Blob at a point
  563.  * is the sum of the gradients of the individual polynomials at that point.
  564.  * The individual polynomials in this case are di(Ri) i = 0,...,num .
  565.  *
  566.  * To find the gradient of the contributing polynomials
  567.  * take di(Ri) and substitute U = Ri^2;
  568.  *
  569.  *      di(U)    = c4i * U^2  +  c2i * U  + c0i
  570.  *
  571.  *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
  572.  *
  573.  *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 
  574.  *
  575.  *        Ux     = 2 * (x-xi) 
  576.  *
  577.  * Rearranging we get
  578.  *  
  579.  *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
  580.  *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
  581.  *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
  582.  * 
  583.  * where
  584.  *         xdi       =   x0-xi
  585.  *         ydi       =   y0-yi
  586.  *         zdi       =   y0-yi
  587.  *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi   
  588.  *       (xi,yi,zi)  is  the center of Metaball i
  589.  *       c4i,c2i,c0i are the coefficients of Metaball i's density function
  590.  */
  591. int
  592. BlobNormal(blob, pos, nrm, gnrm)
  593. Blob *blob;
  594. Vector *pos, *nrm, *gnrm;
  595. {
  596.     register int i;
  597.  
  598.     /*  
  599.      * Initialize normals to zero 
  600.      */
  601.     nrm->x = nrm->y = nrm->z = 0.0;
  602.     /*
  603.      * Loop through Metaballs. If the point is within a Metaball's
  604.      * Sphere of influence, calculate the gradient and add it to the
  605.      * normals
  606.      */
  607.     for(i=0;i < blob->num; i++)
  608.     {
  609.         register MetaVector *sl;
  610.         register Float dist,xd,yd,zd;
  611.  
  612.         sl = &(blob->list[i]);
  613.         xd = pos->x - sl->x;
  614.         yd = pos->y - sl->y;
  615.         zd = pos->z - sl->z;
  616.  
  617.         dist  = xd*xd + yd*yd + zd*zd;
  618.         if (dist <= sl->rs )
  619.         {
  620.             register Float temp;
  621.  
  622.             /* temp is negative so normal points out of blob */
  623.             temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
  624.             nrm->x += xd * temp;
  625.             nrm->y += yd * temp;
  626.             nrm->z += zd * temp;
  627.  
  628.         }
  629.     }
  630.     (void)VecNormalize(nrm);
  631.     *gnrm = *nrm;
  632.     return FALSE;
  633. }
  634.  
  635.  
  636. /*****************************************************************************
  637.  * Calculate the extent of the Blob
  638.  */
  639. void
  640. BlobBounds(blob, bounds)
  641. Blob *blob;
  642. Float bounds[2][3];
  643. {
  644.     Float r;
  645.     Float minx,miny,minz,maxx,maxy,maxz;
  646.     int i;
  647.  
  648.     /*
  649.      * Loop through Metaballs to find the minimun and maximum in each
  650.      * direction.
  651.      */
  652.     for( i=0;  i< blob->num;  i++)
  653.     {
  654.         register Float d;
  655.  
  656.         r = sqrt(blob->list[i].rs);
  657.         if (i==0)
  658.         {
  659.             minx = blob->list[i].x - r;
  660.             miny = blob->list[i].y - r;
  661.             minz = blob->list[i].z - r;
  662.             maxx = blob->list[i].x + r;
  663.             maxy = blob->list[i].y + r;
  664.             maxz = blob->list[i].z + r;
  665.         }
  666.         else
  667.         {
  668.             d = blob->list[i].x - r; 
  669.             if (d<minx) minx = d;
  670.             d = blob->list[i].y - r; 
  671.             if (d<miny) miny = d;
  672.             d = blob->list[i].z - r; 
  673.             if (d<minz) minz = d;
  674.             d = blob->list[i].x + r; 
  675.             if (d>maxx) maxx = d;
  676.             d = blob->list[i].y + r; 
  677.             if (d>maxy) maxy = d;
  678.             d = blob->list[i].z + r; 
  679.             if (d>maxz) maxz = d;
  680.         }
  681.     }
  682.     bounds[LOW][X]  = minx;
  683.     bounds[HIGH][X] = maxx;
  684.     bounds[LOW][Y]  = miny;
  685.     bounds[HIGH][Y] = maxy;
  686.     bounds[LOW][Z]  = minz;
  687.     bounds[HIGH][Z] = maxz;
  688. }
  689.  
  690. char *
  691. BlobName()
  692. {
  693.     return blobName;
  694. }
  695.  
  696. void
  697. BlobStats(tests, hits)
  698. unsigned long *tests, *hits;
  699. {
  700.     *tests = BlobTests;
  701.     *hits = BlobHits;
  702. }
  703.  
  704. void
  705. BlobMethodRegister(meth)
  706. UserMethodType meth;
  707. {
  708.     if (iBlobMethods)
  709.         iBlobMethods->user = meth;
  710. }
  711.